Chemical Relationship among Genetically Authenticated Medicinal Species of Genus Angelica

The genus Angelica comprises various species utilized for diverse medicinal purposes, with differences attributed to the varying levels or types of inherent chemical components in each species. This study employed DNA barcode analysis and HPLC analysis to genetically authenticate and chemically classify eight medicinal Angelica species (n = 106) as well as two non-medicinal species (n = 14) that have been misused. Nucleotide sequence analysis of the nuclear internal transcribed spacer (ITS) region revealed differences ranging from 11 to 117 bp, while psbA-trnH showed variances of 3 to 95 bp, respectively. Phylogenetic analysis grouped all samples except Angelica sinensis into the same cluster, with some counterfeits forming separate clusters. Verification using the NCBI database confirmed the feasibility of species identification. For chemical identification, a robust quantitative HPLC analysis method was developed for 46 marker compounds. Subsequently, two A. reflexa-specific and seven A. biserrata-specific marker compounds were identified, alongside non-specific markers. Moreover, chemometric clustering analysis reflecting differences in chemical content between species revealed that most samples formed distinct clusters according to the plant species. However, some samples formed mixed clusters containing different species. These findings offer crucial insights for the standardization and quality control of medicinal Angelica species.

Previous research endeavors have sought to categorize Angelica species based on their chemical, genetic, or morphological differences.Chromatographic techniques such as TLC, HPLC, and LC/MS have been employed to discern A. sinensis, A. pubescens f. biserrata, A. dahurica, and other related Umbelliferae plants by analyzing coumarins, phthalides, phenolics, and polyacetylene [8].Additionally, seven Angelica species (A.gigas, A. acutiloba, A. tenuissima, A. dahurica, A. koreana, A. polymorpha, and A. decusriva) were authenticated through quantitative analysis of coumarins and micro-morphologies [9].Furthermore, HPLC was utilized to differentiate A. sinensis, A. acutiloba, A. acutiloba var.sugiyamae, and other related Umbelliferae herbs based on their chromatographic fingerprints [10].Lastly, quantitative analysis of coumarins and phenolics using HPLC enabled the chemical differentiation of three Angelica species of Dang-gwi (A. gigas, A. acutiloba, and A. sinensis) [11].
In this study, genetic analysis was employed as a tool to explore the phylogenetic relationships among various Angelica species as well as related plants within the Apiaceae family.These species were genetically classified using a combination of nrDNA internal transcribed spacer (ITS) and external transcribed spacer sequences, cpDNA sequences (rpsl6 intron, rpsl6-tmK, rpl32-trnL, and trnL-trnT), and macro-and micro-morphological characteristics [12].DNA barcoding regions, which included three chloroplast regions (rbcL, matK, and trnH-psbA) and the nuclear ITS region, were utilized to determine phylogenetic relationships among A. sinensis, A. biserrata and A. dahurica [13].Chloroplast genome sequences were used to establish the phylogenetic relationships of 33 Angelica species and 31 other Apioideae species [14].Another study reported the use of 5S-rRNA spacer domains and chemical components (ferulic acid and Z-ligustilide) as genetic and chemical markers, respectively, to compare species differences among A. gigas, A. sinensis, and A. acutiloba [15].However, there were limitations in the study, as the samples of Angelica species used in the chemical analysis were not guaranteed by their exact botanical species, and the chemical relationships among the Angelica samples in the genetic analyses were not confirmed.
As mentioned above, herbal medicines from Angelica species have been utilized for their medicinal purposes in Korean traditional medicine.However, there is controversy in defining the original species of Gang-hwal, of which the botanical origin is Ostericum koreanum Maximowicz in Korean pharmacopeia [1].Therefore, in this study, to find the possible quantitative explanations for the differences in genus Angelica-oriented herbal medicines, we genetically identified eight species of medicinal Angelica genus and two non-medicinal Angelica species, namely, A. polymorpha Maxim.(Korean name: Gunggungi) and Ostericum grossiserratum (Maxim.)Kitag.(=O.koreanum (Maxim.)Kitag., A. grosseserrata Maxim., and A. koreana Maxim.)(Korean name: Singamchae).Those two non-medicinal species were examined to define the original species of Gang-hwal.We used the ITS and psbA-trnH regions for genetic identification.Furthermore, we chemically distinguished these species for chemotaxonomic classification using quantitative HPLC analysis of fortysix marker compounds.Subsequently, we investigated the chemical relationships among medicinal Angelica species as well as non-medicinal species and thereby found feasible alternatives to medicinal species.

DNA Barcode Analysis
To identify the species among 120 samples derived from 8 medicinal species and 2 non-medicinal species of Angelica genus, including 13 samples of A. acutiloba, 10 samples of A. biserrata, 15 samples of A. dahurica, 13 samples of A. decursiva, 11 samples of A. gigas, 24 samples of A. reflexa, 12 samples of A. sinensis, 8 samples of C. tenuissimum, 6 samples of A. polymorpha, and 8 samples of O. grossiserratum, the nucleotide sequences of the ITS and psbA-trnH regions were analyzed.In the ITS region, approximately 688-694 bases of amplified product sequences were examined for each species.No intraspecies variation was observed in either the ITS or psbA-trnH nucleotide sequences.The analysis revealed nucleotide sequence differences ranging from 11 bp to 117 bp depending on the species, with an efficient classification of 10 species (sequence identity matrix range 0.832-0.978,Table S1).For the psbA-trnH region, approximately 307-350 bases of amplified product sequences were analyzed for each species.This region exhibited a 3-95 bp nucleotide sequence difference depending on the species (sequence identity matrix range 0.990-0.738,Table S2).The results of species discrimination in both the ITS and psbA-trnH regions were found to be consistent.The base sequence of the analyzed sample was confirmed for species identification through dual verification using both standard sample data and the NCBI database (Table 1).
. 'Purchased', the samples were purchased from herbal companies in Korea.'Provided (NG)', the samples were provided as the name of 'Nam-Gang-hwal'.'Provided (BG)', the samples were provided as the name of 'Buk-Gang-hwal'.'Collected', the samples were collected from wild habitats.

Phylogenetic Analysis
In the PhyML + SMS (Maximum likelihood-based inference of phylogenetic trees with Smart Model Selection) tree constructed based on concatenated nucleotide sequences of the ITS and psbA-trnH regions (Figure 1), the phylogenetic tree displayed clear separation by species, thereby supporting the accuracy of the identification results based on the two DNA barcode regions.All samples derived from the genus Angelica were clustered together, except for ASI (A. sinense).ASI was closely grouped with L. jeholense and L. tenuissium (=CTE), distinct from other genus Angelica samples.Additionally, one of the non-herbal species, O. grossiserratum (OGR), was closely clustered with other Ostericum species, forming a distinct cluster.
Among the non-specific markers, xanthotoxol (10) and phellopterin (36) exhibited significantly higher levels in the ADA samples compared to both the AAC samples and Plants 2024, 13, 1252 9 of 20 APO samples, respectively.Similarly, the contents of angelol A (23), columbianetin acetate (30), and columbianadin (42) were notably elevated in the ABI samples compared to those in the ADA samples, ARE samples, and ADE samples, respectively.However, no significant quantitative differences were observed for prim-O-glucosyl-cimifugin (3) between the ADA and ARE samples nor for osthol (37) between the ABI and ARE samples.
Nodakenin (5) exhibited significantly higher contents in both the ADE and AGI samples, while umbelliferone (6) displayed elevated levels specifically in the ADE samples and benzoic acid (8) demonstrated increased content in both the ASI and CTE samples compared to the AAC samples.Moreover, the contents of marmesin (12), decursinol (15), decursin (38), and decursinol angelate (39) were notably higher in the AGI samples.Similarly, byakangelicin (17), byakangelicol (28), and imperatorin (34) exhibited elevated levels in the ADA samples, while xanthotoxin (22) was more abundant in the AAC samples.Additionally, bergapten (25) showed increased content in the ABI samples, while ostenol (26) and bisabolangelone (27) displayed higher levels in the ARE samples.Moreover, senkyunolide A (32) exhibited elevated content in the CTE samples and falcarindiol (43) demonstrated significantly higher levels in the OGR samples compared to other samples.
Chlorogenic acid (1) displayed significantly higher levels in the AGI and ARE samples, while ferulic acid (7) exhibited elevated content in the ARE and ASI samples.Senkyunolide I ( 9), senkyunolide H (11), ligustilide (35), and levistilide A (45) demonstrated increased levels in the ASI and CTE samples, while oxypeucedanin hydrate (14) and isoimperatorin (40) showed elevated contents in the ADA samples.Additionally, coniferyl ferulate (31) exhibited significantly higher levels in the APO and CTE samples compared to other samples.In contrast, caffeic acid (2) showed significantly lower contents in the ASI and CTE samples and 3-n-butyl-phthalide (33) demonstrated decreased levels in the AAC samples compared to other samples.

Chemometric Clustering Analysis
In the hierarchical clustering analysis (HCA), the majority of samples formed distinct clusters corresponding to their botanical species, with the occasional insertion of samples from other species.Notably, AAC samples (excluding AAC01, AAC06, and AAC12) and ABI05 clustered together, along with a close positioning of a mixture of samples (AAC12, ABI01, -02, ASI02, -03, -07, and CTE02).The remaining ABI samples were clustered with AAC06, ARE14, and ARE16.ADA samples formed a separate cluster without any interspersions.ADE samples were split into two distinct clusters within their species (ADE02−ADE08 vs. ADE09-ADE13).All AGI samples, along with AAC01 and ARE07, formed a distinct cluster.APO samples, except for APO04, also formed their own cluster.Although three samples (ARE07, ARE14, and ARE16) were positioned adjacent to other samples, ARE samples were divided into two separate clusters, closely associated with ADE samples and APO and ADA samples, respectively.The remaining ASI samples were combined with CTE01 and CTE07, near the cluster of CTE samples.OGR samples formed a separate cluster, distinctly different from other samples (Figure 3).
In the principal component (PC) score plot, most ADA samples and ARE samples exhibited negative PC1 scores and positive PC2 scores.Conversely, ASI and CTE samples displayed positive scores for both PC1 and PC2.The ABI samples, except for two samples, had negative PC2 scores but varied in their PC1 scores between positive and negative values.These sample distributions were distinctly discernible based on PC scores.However, there were compact distributions of AAC, ADE, AGI, APO, OGR, and some ARE samples near 'zero' PC scores.Similar to the clustering in Figure 3, the ARE samples were distributed as a Nam-Gang-hwal group including ARE07 and -20 and a Buk-Gang-hwal group including ARE13 and -14.Although samples were grouped by their respective species, there was notable overlap in their distributions (Figure 4).formed a distinct cluster.APO samples, except for APO04, also formed their own cluster.Although three samples (ARE07, ARE14, and ARE16) were positioned adjacent to other samples, ARE samples were divided into two separate clusters, closely associated with ADE samples and APO and ADA samples, respectively.The remaining ASI samples were combined with CTE01 and CTE07, near the cluster of CTE samples.OGR samples formed a separate cluster, distinctly different from other samples (Figure 3).In the principal component (PC) score plot, most ADA samples and ARE samples exhibited negative PC1 scores and positive PC2 scores.Conversely, ASI and CTE samples displayed positive scores for both PC1 and PC2.The ABI samples, except for two samples, had negative PC2 scores but varied in their PC1 scores between positive and negative values.These sample distributions were distinctly discernible based on PC scores.However, there were compact distributions of AAC, ADE, AGI, APO, OGR, and some ARE samples near 'zero' PC scores.Similar to the clustering in Figure 3, the ARE samples were distributed as a Nam-Gang-hwal group including ARE07 and -20 and a Buk-Gang-hwal group including ARE13 and -14.Although samples were grouped by their respective species, there was notable overlap in their distributions (Figure 4).

Correlation Analysis
Chemical correlations among Angelica species were assessed by computing Pearson's correlation coefficient (r) for both intra-and interspecies comparisons, as outlined in Table 2 and illustrated in Figure 5.Among the Angelica species, AGI, APO, ASI, CTE, and OGR samples exhibited notably high intraspecies coefficients, with mean and median r values > 0.9.Following closely were the ADA and AAC samples, displaying mean r values ranging from 0.6 to 0.8 and median r values ranging from 0.8 to 0.9.However, distinct outliers were observed among the samples such as AAC01, AAC06, ASI03, ASI07, and CTE02 in terms of intraspecies coefficients.In contrast, ADE and ARE samples showed lower mean r values of 0.4−0.5 and median r values of 0.3−0.4,respectively.The ABI samples exhibited the lowest mean and median r values of coefficients, both below 0.2 (Figure S3).

Correlation Analysis
Chemical correlations among Angelica species were assessed by computing Pearson's correlation coefficient (r) for both intra-and interspecies comparisons, as outlined in Table 2 and illustrated in Figure 5.Among the Angelica species, AGI, APO, ASI, CTE, and OGR samples exhibited notably high intraspecies coefficients, with mean and median r values > 0.9.Following closely were the ADA and AAC samples, displaying mean r values ranging from 0.6 to 0.8 and median r values ranging from 0.8 to 0.9.However, distinct outliers were observed among the samples such as AAC01, AAC06, ASI03, ASI07, and CTE02 in terms of intraspecies coefficients.In contrast, ADE and ARE samples showed lower mean r values of 0.4-0.5 and median r values of 0.3-0.4,respectively.The ABI samples exhibited the lowest mean and median r values of coefficients, both below 0.2 (Figure S3).
Interspecies correlations presented a wider range of values compared to intraspecies correlations.AAC samples displayed relatively higher mean and median coefficients with ASI and CTE samples (both r > 0.5), while coefficients with other samples were generally below 0.2 (except for ARE samples with mean r values > 0.2) and even showed negative values with ADA samples.ABI samples showed mean coefficients below 0.2, with most median values being negative.ADA samples exhibited mean and median r values ranging from 0.2 to 0.4 with APO and ARE samples but negative values with samples from other species.ADE samples displayed relatively higher mean coefficients with APO samples (r > 0.4) and OGR samples (r > 0.5) and lower mean values with ARE samples (r > 0.1).AGI samples showed coefficients mostly close to zero, with both positive and negative values, when compared to other samples.APO samples showed higher coefficients with OGR samples (mean r value > 0.8 and median r value > 0.7), followed by higher values with ARE samples (mean and median r values > 0.4).However, the coefficients of CTE samples with OGR samples were negative in both mean and median values (Figure S4).

Discussion
Recently, entire cp genome analysis has been widely used for gene-based species discrimination or identification research.Those versatile characteristics make the cp region specifically applicable to plants, and its genetically informative feature is quite an attractive approach to botanical identification.However, the cp genome analysis has a limitation on the samples being processed via multiple steps.Therefore, DNA barcoding and phylogenetic analysis were used to identify the specific botanical species among samples of the Angelica genus.Previous studies highlighted the effectiveness of ITS as a robust tool, but analysis of the psbA-trnH region was also conducted in this study to increase result accuracy [22,23].Despite some reported issues with indels in the psbA-trnH region Pearson's correlation coefficient (r) analysis further confirmed the stronger associations of Nam-Gang-hwal samples with ADA samples (mean r value = 0.44) and Buk-Gang-hwal samples with AAC (mean r value = 0.34) and ADE samples (mean r value = 0.34), compared to the correlations of Buk-Gang-hwal samples with ADA samples (mean r value = 0.09) and Nam-Gang-hwal samples with AAC (mean r value = 0.13) and ADE samples (mean r values < 0), respectively.Meanwhile, APO and OGR samples displayed closer correlations with Buk-Gang-hwal samples (both mean r value > 0.5) than with Nam-Gang-hwal samples.

Discussion
Recently, entire cp genome analysis has been widely used for gene-based species discrimination or identification research.Those versatile characteristics make the cp region specifically applicable to plants, and its genetically informative feature is quite an attractive approach to botanical identification.However, the cp genome analysis has a limitation on the samples being processed via multiple steps.Therefore, DNA barcoding and phylogenetic analysis were used to identify the specific botanical species among samples of the Angelica genus.Previous studies highlighted the effectiveness of ITS as a robust tool, but analysis of the psbA-trnH region was also conducted in this study to increase result accuracy [22,23].Despite some reported issues with indels in the psbA-trnH region due to its non-coding region, it remains an efficient tool for analyzing herbal medicines due to its relatively short length and abundant variation compared to other cpDNA barcode regions such as rbcL and MatK [24].The results of this study also revealed that while the psbA-trnH sequence length is approximately half that of ITS, the variation in sequence between species is similar to that of ITS.
The taxonomic classification of A. reflexa (Korean name: Gang-hwal), of which the Korean botanical name is the same as the Korean herbal name, has been quite complex.It was initially classified as A. koreana by Maximowicz (1886) [25] and later transferred to the genus Ostericum (O.koreanum) by Kitagawa (1936) [26] due to its external morphological similarity.Subsequently, Kitagawa (1971) [27] recognized the taxonomic identity between A. koreana (=O.koreanum) and O. grossiserratum, treating A. koreana (=O.koreanum) as a synonym of O. grossiserratum.This is the reason why O. koreanum still remains a synonym of O. grossiserratum.However, molecular phylogenetic study indicated that A. koreana (=O.koreanum) may be independent from O. grossiserratum [28].Furthermore, both external morphological examination and molecular phylogeny using nuclear DNA ITS sequences have shown that commercial medicinal plants cultivated as Gang-hwal are neither A. koreana nor O. grossiserratum [29].After careful observation of morphological and anatomical characters and examination of relevant specimens, Lee et al. (2013) [5] ultimately proposed this as a new species of Angelica, named A. reflexa.Some researchers still consider Gang-hwal to be A. genuflexa due to strong morphological similarity and regard A. reflexa as a synonym of A. genuflexa [29].However, careful observation [5] and our DNA barcode analysis support the difference between A. reflexa and A. genuflexa.As shown in Figure 1, all A. reflexa (ARE) samples were located within the 'Angelica' clade, not within the 'Ostericum' clade, while O. grossiserratum was clearly classified as genus Ostericum and separated from A. genuflexa.According to the results of phylogenetic analysis, A. genuflexa was closely grouped with A. biserrata (ABI) and A. polymorpha (APO) rather than A. reflexa.Moreover, the herbal samples of Gang-hwal acquired for this study were clearly identified as A. reflexa, not O. koreanum or O. grossiserratum.
Interesting findings were observed for the three Angelica species commonly used as Dang-gwi (Figure 1).Despite A. gigas (AGI) and A. acutiloba (AAC) being grouped within the same Angelica clade, they are distinctly separated into different clusters (AGI clustered with ADE).Additionally, A. sinensis (ASI) was clustered with L. sinense and L. jeholense, positioned within the Ligusticum clade including CTE samples rather than the Angelica clade.
The roots of A. reflexa are utilized in Korean herbal markets to produce two distinct types of Gang-hwal, known as 'Nam-Gang-hwal' through seed-propagation and 'Buk-Gang-hwal' via root propagation [40,41].These differing cultivation methods yield variations in the chemical compositions between Nam-Gang-hwal and Buk-Gang-hwal [35].In this study, the ARE samples were also categorized into the Nam-Gang-hwal group and the Buk-Gang-hwal group.These differences in sample types led to the division of ARE samples into separate groups in both the HCA dendrogram and the PC score plot.Additionally, Nam-Gang-hwal samples exhibited a closer chemical relationship with the ADA samples, while Buk-Gang-hwal samples showed a closer relationship with the ADE and AAC samples, forming clusters of closely related samples [42].
The chemical relationships among three Angelica species of Dang-gwi were depicted differently in the dendrogram and Pearson's coefficients.In the dendrogram, the AAC samples and ASI samples were distinctly separated into individual clusters, while the AGI samples formed a secondary cluster, situated apart from the clusters of AAC and ASI samples.Pearson's coefficients also highlighted the chemical heterogeneity of AGI samples from AAC and ASI samples, with mean r values < 0.2.However, a higher r-value of AAC-ASI (mean r value > 0.5) indicated a closer chemical relationship between AAC and ASI samples compared to the combination with AGI samples.Previous studies have reported chemical heterogeneity among three Angelica species of Dang-gwi, but in contrast to these findings, A. gigas and A. acutiloba exhibited a closer relationship than with A. sinensis samples [10,14].
Overall, the chemical relationship between ASI and CTE samples appears to reflect their phylogenetic relevance, despite their difference in therapeutic activities.The ARE samples (i.e., Nam-Gang-hwals) exhibited chemical and phylogenetic relevance with ADA samples, while ARE samples (i.e., Buk-Gang-hwals) showed relevance with AAC and ADE samples.The APO samples, classified as non-medicinal species, displayed partial chemical and genetic relevance with ABI, ADA, and ARE samples and a stronger correlation with OGR samples.Although the OGR samples, another non-medicinal species, exhibited the lowest genetic relationship with other species, their chemical relationship with other species varied, either parallel or opposite to the genetic results, depending on the statistical tools used.
This study has several limitations: (1) the unequal distribution of samples among species; (2) the limited representativeness of the selected 46 marker compounds to whole chemical characteristics of all species samples; (3) inconsistencies in the chemical relationships among the samples in the chemometric analyses; and (4) the lack of comparison between ARE samples and their therapeutic analogous Notopterygium species.Despite these limitations, this study represents the first attempt to genetically authenticate and chemically classify medicinal Angelica species, as well as non-medicinal species.The quantitative explanation of the differences among the Angelica species, including the alternative use of non-medicinal species, would be further supported by pharmacological study.

Preparation of Genomic DNA
The genomic DNA was extracted from the samples following the instructions provided in the NucleoSpin ® Plant II kit manual (Macherey-Nagel, Dueren, Germany).This process involved the use of a PL1 lysis buffer during a lysis step that lasted at least 2 h.For certain samples, an additional step was incorporated, which involved the use of 10% cetyltrimethyl ammonium bromide (CTAB) and 0.7M NaCl.This extra step was employed to remove phenolic compounds and polysaccharides after the extraction of DNA using the kit.

Preparation of Genomic DNA
The genomic DNA was extracted from the samples following the instructions provided in the NucleoSpin ® Plant II kit manual (Macherey-Nagel, Dueren, Germany).This process involved the use of a PL1 lysis buffer during a lysis step that lasted at least 2 h.For certain samples, an additional step was incorporated, which involved the use of 10% cetyltrimethyl ammonium bromide (CTAB) and 0.7M NaCl.This extra step was employed to remove phenolic compounds and polysaccharides after the extraction of DNA using the kit.

Determination of DNA Sequences of PCR Product
The PCR product, separated from the agarose gel using the MagListo™ 5M PCR/Gel Purification Kit (Bioneer, Daejeon, Republic of Korea), was cloned using the TOPcloner™ TA Kit (Enzynomics, Daejeon, Republic of Korea).The DNA sequences of the cloned PCR product were then determined through analysis performed by Bioneer (Daejeon, Republic of Korea).

Analysis of DNA Sequences and Preparation of Phylogenetic Tree
DNA sequences were analyzed using ClustalW multiple sequence alignment (Bioedit, v7.7.1) and confirmed with multiple sequence alignment in MAFFT (MAFFT, v7) [45].To verify the polymorphisms represented by IUPAC symbols in the sequence data, all sequences were generated at least twice.The chromatograms of nucleotide sequences, provided by the Bioneer sequencing service, were compared.Evolutionary analyses were conducted in MEGA X (ver.10.0.5).Phylogenetic analysis of two concatenated DNA barcode regions (ITS and psbA-trnH) was constructed using the PhyML + SMS/OneClick method, which showed a workflow of MAFFT, BMGE, and PhyML + SMS (maximum likelihood-based inference of phylogenetic trees with Smart Model Selection) [46].All analyzed sequences were compared with NCBI GenBank using BLAST [47].Newly determined nucleotide sequences were deposited in NCBI GenBank.Two other subfamilies of the Apiaceae, Eryngium regnellii (Saniculoideae) and Centella asiatica (Mackinlayoideae), were used as outgroups.NCBI data used in the phylogenetic tree analysis were represented with the accession number and scientific names listed in the NCBI database.

Analytical Sample Preparation
Before use, all samples (each triplicate) were thoroughly dried and then ground to a powder, which was homogenized through a 500 µm testing sieve (Chunggyesanggong-sa; Gunpo, Republic of Korea).A precise weight of the powder (500 mg) was extracted with 5 mL of methanol for 30 min using an ultrasonic extractor (Power Sonic 520; Hwashin Tech, Daegu, Republic of Korea).The extract was then centrifuged at 10,000 rpm for 10 min.The supernatant was transferred to a 1.5-mL microtube and gently dried using a nitrogenblowing concentrator (MGS2200; Eyela, Tokyo, Japan).The residue was re-dissolved in HPLC-grade methanol at a concentration of 10,000 µg/mL and filtered through a 0.2 µm syringe filter (BioFact, Daejeon, Republic of Korea) before HPLC injection.

HPLC Analytical Conditions
The quantitative analysis of the marker compounds was performed using an Agilent 1200 liquid chromatography system (Agilent Technologies, Palo Alto, CA, USA), equipped with an autosampler, degasser, quaternary solvent pump, and diode array detector.The data acquired were processed using Chemstation software (Rev.B. 04.03.; Agilent Technologies Inc., USA).A Capcell Pak Mg II C 18 column (4.6 mm × 250 mm, 5 µm; Shiseido, Tokyo, Japan) was used to separate forty-six marker compounds at 35 • C with a flow rate of 1 mL/min and an injection volume of 10 µL.The mobile phases were pumped via gradient elution by mixing water containing 0.1% TFA (solvent A) and acetonitrile (solvent B).The percentages of mobile phase (solvent B) with equivalent retention times were as follows: 15% for 0-2 min, 15-50% for 2-30 min, 50-50% for 30-32 min, 50-75% for 32-55 min, and 75% for 55-58 min, and then re-equilibrated to 15% until the end of the analysis.The detection wavelength of the diode-array detector was set at 230, 250, 270, 280, and 325 nm.

Validation of the HPLC Method
The stock solution was prepared by dissolving each marker compound in methanol at a concentration of 1000 µg/mL, and the working solution for the calibration curve was generated by serial dilution of the stock solution to seven different concentrations.The correlation coefficients (r 2 ) were determined to assess the linearity of the calibration curve.The limit of detection and the limit of quantification were established as signal-to-noise (S/N) ratios of 3 and 10, respectively.
Precision, indicative of the repeatability of the analytical method, was assessed by analyzing low and high concentrations of the stock solutions three times within one day (intraday precision) and over three consecutive days (interday precision).Precision was expressed as relative standard deviations (RSDs): RSD (%) = (standard deviation/ mean) × 100.
The accuracy of the analytical method, represented by recovery, was evaluated by adding low and high concentrations of the marker compounds to the sample solutions.The equation for calculating recovery was as follows: Recovery (%) = [(detected concentrationinitial concentration)/spiked concentration] × 100.

Table 1 .
Species identification of the samples based on DNA barcode analysis of the ITS and psbA-trnH regions.

Table 2 .
Summary of the Pearson's correlation coefficients of the samples.